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Water  management  in  polymer  electrolyte  membrane  (PEM)  fuel  cells  is  important  for  fuel  cell  perfor¬ 
mance  and  durability.  Numerical  simulations  using  the  lattice  Boltzmann  method  (LBM)  are  developed 
to  elucidate  the  dynamic  behavior  of  condensed  water  and  gas  flows  in  a  PEM  fuel  cell  gas  channel.  A 
scheme  for  two-phase  flow  with  large  density  differences  was  applied  to  establish  the  optimum  gas  chan¬ 
nel  design  for  different  gas  channel  heights,  droplet  initial  positions,  droplet  volume  and  air  flow  velocity 
for  both  hydrophobic  and  hydrophilic  gas  channels.  The  discussion  of  optimum  channel  height  and  drain 
performance  was  made  using  two  factors  “pumping  efficiency”  and  “drainage  speed”.  It  is  shown  that 
deeper  channels  give  better  draining  efficiency  than  shallower  channels,  but  the  efficiency  dramatically 
decreases  when  the  droplet  touches  corners  or  the  top  of  gas  channel’s  walls.  As  the  droplet  velocity,  i.e. 
the  drainage  flow  rate  becomes  higher  and  the  drainage  efficiency  becomes  less  dependent  on  droplet 
locations  with  shallower  channels,  shallower  channels  are  better  than  deeper  channels.  Introducing  a 
new  dimensionless  parameter,  “pumping  efficiency”,  the  investigation  discusses  the  effect  of  the  various 
parameters  on  the  drainage  performance  of  a  PEM  fuel  cell  gas  channel. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Phenomena  involved  in  PEM  fuel  cell  operation  are  com¬ 
plex;  specifically,  they  involve  heat  transfer,  species  and  charge 
transport,  multiphase  flows,  and  electrochemical  reactions.  Basic 
research  into  these  phenomena  is  critically  important  to  over¬ 
come  two  major  barriers  to  PEM  fuel  cell  use,  durability  and 
cost.  This  paper  investigates  water  management  which  is  essential 
for  improving  the  performance  of  polymer  electrolyte  membrane 
(PEM)  fuel  cells.  The  membrane  of  PEM  fuel  cells  has  to  be  fully 
hydrated  to  maintain  high  proton  conductivity,  and  at  the  same 
time  excess  water  condenses  in  the  gas  diffusion  layers  (GDLs)  or 
gas  flow  channels  (GFCs)  and  prevents  the  supply  of  reactants  to 
the  electrodes  under  high  current  density  conditions.  Phenomena 
related  to  this  are  generally  referred  to  as  “flooding”  and  may  be 
a  cause  of  durability  and  performance  reductions  due  to  reactant 
starvation,  and  the  GDL  generally  uses  hydrophobic  materials  to 
facilitate  liquid  water  drainage,  like  in  the  investigation  of  the  LBM 
simulations  reported  here.  At  the  cathode  GDL/GFC  interface,  oxy¬ 
gen  transports  towards  the  electrode  where  it  reacts  with  protons 
and  electrons  to  produce  water,  which  eventually  enters  the  cell 


*  Corresponding  author.  Tel.:  +81  11  706  6333;  fax:  +81  11  706  7889. 
E-mail  address:  yasser@eng.hokudai.ac.jp  (Y.  Ben  Salah). 


channels.  The  interfacial  resistance  to  the  reactant  transport  will  be 
significantly  increased  by  the  presence  of  liquid  water  here.  Opti¬ 
cal  visualization  has  shown  that  liquid  water  is  present  as  droplets 
on  the  GDL  surface,  and  is  removed  by  the  gas  flow  and/or  attach¬ 
ing  to  the  channel  walls  [1].  Studies  have  been  conducted  on  the 
liquid  water  behavior  in  channels  and  optimization  of  gas  channel 
design.  Chen  et  al.  [2]  conducted  the  analysis  of  droplet  instabil¬ 
ity  and  detachment  and  indicated  that  the  static  contact  angle 
{0S)  and  contact  angle  hysteresis  (the  difference  between  advanc¬ 
ing  and  receding  contact  angles,  i.e.  0A  -0R),  are  both  important 
parameters  in  determining  the  force  required  to  move  a  droplet 
across  a  surface.  Instability  diagrams  were  developed  to  explore 
the  operating  conditions  under  which  droplets  become  unstable, 
as  unstable  droplet  conditions  are  desirable  to  operate  the  fuel  cell 
under  conditions  allowing  the  instantaneous  removal  of  droplets 
from  the  GDL/GFC  interface  so  as  to  prevent  blockage  of  pathways 
for  oxygen  transport  to  the  three-phase  reaction  sites.  Following 
Chen’s  work  [2],  Kumbur  et  al.  [3]  proposed  a  similar  analytical 
model  to  elucidate  the  effects  of  channel  geometry  and  GDL  surface 
properties  on  water  droplet  instability.  Hao  et  al.  [4]  used  the  mul¬ 
tiphase  free-energy  lattice  Boltzmann  method  to  study  the  effect 
of  gas  flow  velocity  and  GDL  wettability  on  water  droplet  dynamic 
behavior.  Two-dimensional  simulation  employing  the  volume  of 
fluid  (VOF)  method  were  performed  to  investigate  the  dynamic 
behavior  of  a  water  droplet  subjected  to  air  flow  in  the  bulk  of 
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Nomenclature 

c  characteristic  particle  speed  (m  s-1 ) 

Cj  restricted  velocities  of  particle  ensembles  (m  s-1 ) 

fi  particle  velocity  distribution  functions  for  the  cal¬ 

culation  of  an  order  parameter 
g  gravitational  acceleration  (ms-2) 

gi  particle  velocity  distribution  functions  for  the  cal¬ 

culation  of  a  predicted  velocity 
H  vertical  length  of  simulation  domain  (m) 

I  cell  current  density  (A  m-2 ) 

L  characteristic  length  (m) 

m  droplet  mass  (kg) 

p  pressure  (Pa) 

Q  gas  flow  rate  (m3  s-1 ) 

Sh  Strouhal  number 

t  time  (s) 

t0  characteristic  time  scale  (s) 

At  time  step  during  which  the  particles  travel  the  lat¬ 
tice  spacing  (s) 

u  current  velocity  (ms-1) 

u*  predicted  velocity  (m  s-1 ) 

U  characteristic  flow  speed  (m  s-1 ) 

Udr  droplet  velocity  (ms-1) 

x,y,z  position  coordinates  (m) 

Ax  spacing  of  the  cubic  lattice  (m) 

Kf  constant  determining  the  width  of  the  interface  of 

two  phases 

Kg  constant  determining  the  strength  of  the  surface 
tension 

0  contact  angle 

fi  viscosity  (Pas) 

§  coordinate  perpendicular  to  the  interface  (m) 

p  density  (kg  m-3) 

p0  reference  density  ( kg  m-3 ) 

a  interface  tension  (N  m-1 ) 

Tf,Tg  dimensionless  single  relaxation  time 

rj  pumping  efficiency 

0  order  parameter 

0o  reference  order  parameter 

Superscripts/subscripts 
A  advancing 

dr  droplet 

eq  equilibrium  state 

ext  exterior 

in  inflow 

int  interior 

G  gas 

I  liquid 

R  receding 

5  solid 

s  static 

sd  difference  between  with  and  without  droplet 

oi,/3  Cartesian  coordinates 


the  gas  channel  [5]  and  to  study  the  detachment  of  liquid  droplets 
from  the  surfaces  of  porous  materials  used  in  (PEM)  fuel  cells, 
under  the  influence  of  cross-flowing  air  [6].  The  effects  of  gas 
flow  velocity  and  surface  wettability  on  the  two-phase  flow  pat¬ 
terns  in  flow  channels  were  investigated  most  recently  by  Ding 
et  al.  [7]  using  the  volume  of  fluid  (VOF)  method.  The  VOF  method 
could  include  the  effect  of  dynamic  contact  angle  changes,  which 
is  an  important  parameter  in  the  droplet  dynamics  in  the  present 


model.  In  order  to  obtain  the  dynamic  change  of  the  contact  angle, 
a  complicated  numerical  scheme  must  be  used  to  track  interface 
changes  continuously  [6],  and  experimental  correlations  for  the 
advancing  and  receding  contact  angles  with  the  gas  velocity  were 
employed  in  the  VOF  method  [8].  The  lattice  Boltzmann  method 
(LBM)  is  a  powerful  technique  for  simulating  transport  and  fluid 
flows  involving  interfacial  dynamics  and  complex  geometries.  In 
particular,  due  to  the  kinetic  nature  and  absence  of  a  need  to  track 
the  phase  interfaces,  the  LBM  has  been  found  very  effective  to  sim¬ 
ulate  two-phase  flow  in  the  gas  channels  [9,10].  The  LBM  could 
also  estimate  the  relation  between  the  dynamic  contact  angles  and 
the  droplet  motion  using  the  static  nature  of  wettability  without 
further  experimental  correlations  (as  it  will  be  discussed  later  in 
Section  3.1). 

In  this  paper,  numerical  simulation  using  the  LBM  has  been 
developed  with  for  conditions  with  large  density  differences,  to 
understand  the  dynamic  behavior  of  liquid  water  in  gas  chan¬ 
nels  and  the  effect  of  different  parameters  on  the  draining 
performance.  The  effect  of  droplet  position,  surface  wall  wetta¬ 
bility,  and  channel  height  under  a  constant  flow  rate  are  also 
discussed. 


2.  Simulation  method 

The  LBM  simulates  mass  and  heat  transport  phenomena 
by  tracking  movements  of  particle  ensembles  with  velocities 
restricted  to  a  finite  set  of  vectors.  The  particle  population  is 
expressed  by  distribution  functions,  and  the  time  evolution  of 
the  distribution  functions  is  calculated  by  the  simple  law  of 
collision  and  transition,  ensuring  that  the  LBM  obeys  the  continu¬ 
ity  equation  and  the  Navier-Stokes  equations  for  incompressible 
fluids.  Additionally,  introducing  interaction  among  the  parti¬ 
cle  ensembles  in  the  equations  makes  it  possible  to  simulate 
multi-phase  flow.  Because  of  the  simplicity  of  the  algorithm, 
the  LBM  has  the  following  advantages:  flexibility  for  complex 
boundary  geometries,  simplicity  of  parallel  computing,  and  accu¬ 
racy  in  mass  conservation.  In  multi-phase  flows,  no  tracking 
of  interfaces  is  necessary  and  the  clearly  distinguishable  inter¬ 
faces  can  be  maintained  without  additional  assumptions.  To 
simulate  the  two-phase  flow  in  the  3-dimensional  gas  chan¬ 
nel  of  a  PEM  fuel  cell,  the  extended  LBM  proposed  by  Inamuro 
et  al.  [11]  was  applied.  Two-phase  flows  with  large  density  dif¬ 
ferences,  density  ratios  up  to  1000,  can  be  calculated  by  this 
method  [10]. 

In  the  model,  the  non-dimensional  variables  defined  by  a  char¬ 
acteristic  length  L,  a  characteristic  particle  speed  c,  a  characteristic 
time  scale  t0  =L/U,  where  U  is  a  characteristic  flow  speed,  a  refer¬ 
ence  order  parameter  0O,  and  a  reference  density  p0  are  also  used 
[11];  “non-dimensional”  terms  are  represented  by  a  circumflex. 
This  paper  uses  a  three-dimensional  15  velocities  model  (3D15V 
model)  and  the  velocities  of  particle  ensembles  are  restricted  to 
the  following  vectors  cz  (i  =  1,  2, . . .,  15)  in  the  3-dimensional  case 
as  shown  in  Fig.  1  [12]: 


[Cl ,  C2,  C3,  C4,  C5,  C6,  C7,  Cg,  Cg,  C\q,  C\\ ,  C\2,  C 13,  C 14,  C15] 


"0  1  0  0  -1  0  0  1  -1  1  1  -1  1  -1-1" 
0  0  1  0  0  -1  0  1  1  -1  1  -1  -1  1  -1 

0  0  0  1  0  0  -1  1  1  1  -1  -1  -1  -1  1 


(1) 


Two  particle  velocity  distribution  functions,/)  and  git  are  used. 
The  fi  function  is  used  for  the  calculation  of  an  order  param¬ 
eter  0  which  distinguishes  two  phases:  0  <  0G  corresponding 
to  the  gas  phase,  0  >  0l  the  liquid  phase,  and  0G  <  0  <  0l  the 
condition  at  the  interface  between  liquid  and  gas  phases.  The 
gi  function  is  used  for  the  calculation  of  a  predicted  velocity 
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C14  C\\ 


r  3  9  3 

gjq  =  Ei  1  +3 CiaUa  -  -UaUa  +  +  j  (% 

r 

\  (  dU/3  dua  \  ^  ^ 

)  Ax  ( Wa  +  Wp  )  c,aC,p_ 

+  Ei  p  Gap(p)Cia cifj 

<1 

KJ 

(7) 

where 

E2  =  E3  =  ...  =  E7  =  1, 

Es  =  Eq  =  •  •  •  =  £15  = 

1 

72 

Hi  =  1, 

H2  =  H3  =  •  •  ■  =  H15  =  0, 

E,  =  3E,-  (i  =  2,3, , 

15) 

(7) 

Fig.  1.  Lattice  structure  of  the  three-dimensional  3D15V  LBM. 


for  the  two-phase  fluid  without  a  pressure  gradient.  The  evo¬ 
lution  of  the  particle  distribution  functions  ft  and  g,  with  the 
velocity  cz  at  point  x  and  time  t  are  computed  by  the  following 
equations. 


ft{k  +  CjAt,  t  +  At)  =ft{x,  t)  - 


Tf 


(2) 


gi{x  +  CjAt,  t  +  At)  =  gfx,  t)  -  -j-  [&•(*,  t)  -  gfq(x,  t)] 


+  3TjCja 


a 


/X 


[ft  ,  A 1 
«  y  j 


Ax 


Pc 


-  3£jCjy  1  -  ^  gAx 


(3) 


Here,  f^q  and  gfq  are  the  equilibrium  distribution  functions,  Ez 
is  the  associated  weight  coefficients  presented  below,  z y  and  rg 
are  dimensionless  single  relaxation  times,  Ax  is  the  spacing  of 
the  cubic  lattice,  At  is  the  time  step  during  which  the  parti¬ 
cles  travel  the  distance  of  the  lattice  spacing,  p  is  the  density, 
jl  is  the  viscosity,  u  is  the  current  velocity,  and  g  is  the  gravita¬ 
tional  acceleration.  The  third  and  last  terms  on  the  right  hand  side 
of  Eq.  (3)  represent  the  effects  of  viscous  stress  and  gravitation, 
respectively. 

The  order  parameter  0  distinguishing  the  two  phases  and  the 
predicted  velocity  it  of  the  multi-component  fluids  are  defined  in 
terms  of  the  two  particle  velocity  distribution  functions  as  follows: 


15 

Eft 

(4) 

1=1 

15 

= 

1=1 

(5) 

The  equilibrium  distribution  functions/^9  andg^  in  Eqs.  (2)  and 
(3)  are  given  by 


fieq=Hij>  +  Fi 


A  ^  d2(j)  Kf  f  d(j) 

e(  as 


T  3Ej0Cjaua  +  EiKfGap((j))CiaCip 


(6) 


with  a,  p  =  x,  y,  z  the  subscripts  a  and  /3  represent  the  Cartesian 
coordinates  and  the  summation  convention  is  used.  In  the  above 
equations,  kj  is  a  constant  parameter  determining  the  width  of  the 
interface  between  two  phases,  Kg  is  a  constant  parameter  determin¬ 
ing  the  strength  of  the  surface  tension,  and  the  parameters  p0  and 
Gap  are  explained  in  Ref.  [13].  The  interfacial  tension  a  is  obtained 
by  the  following  equation: 


Here,  |  is  the  coordinate  perpendicular  to  the  interface. 

Because  the  predicted  velocity  it  given  by  Eq.  ( 5 )  does  not  satisfy 
the  continuity  equation  (V  •  u*  =  0),  a  correction  of  it  is  required. 
The  current  velocity  u  which  satisfies  the  continuity  equation  can 
be  obtained  with  the  following  equations: 


Sh 


u-u 

At 


Vp 

P 


(9) 


=  Sh 


V-u* 

At 


(10) 


Here,  Sh  =  U/c  is  the  Strouhal  number  and  p  is  the  pressure  of  the 
two-phase  fluid;  note  that  this  definition  leads  to  the  following 
relationships,  At  =  Sh  Ax,  which  is  represented  by  A t=  A x/c  and 
means  that  the  particles  travel  across  the  lattice  space  Ax  during 
time  step  At.  This  paper  solved  Eq.  (10)  using  the  Successive  Over 
Relaxation  (SOR)  method.  Details  of  this  model  are  described  in  a 
previous  paper  [13]. 

The  effect  of  wettability  is  introduced  by  assuming  the  density  of 
the  solid  wall  as  proposed  by  the  scheme  of  Seta  and  Takahashi  [14]. 
Since  the  intermolecular  force  is  expressed  in  terms  of  the  density 
of  the  fluid  in  the  LBM,  giving  the  density  of  solid  wall  corresponds 
to  giving  the  intermolecular  force  between  liquid  and  solid  wall. 
It  has  been  confirmed  that  this  scheme  can  simulate  the  effect  of 
wettability  both  on  a  flat  surface  as  well  as  at  corners  inside  a  gas 
channel  [15]. 


3.  Results  and  discussion 

3.1.  Basic  characteristics  of  droplets  with  the  LBM 

The  lattice  Boltzmann  method  (LBM)  for  two-phase  flows  with 
large  density  differences  has  been  applied  to  the  simulation  of  liq¬ 
uid  water  and  air  flow  in  a  PEM  fuel  cell  [13].  However,  there  were 
problems  with  the  applicability  of  the  simulation  results,  e.g.  there 
remained  the  issue  of  non  conservation  of  the  mass  of  liquid  water. 
Improvements  to  the  calculation  process,  the  formulations  for  the 
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Fig.  2.  Droplet  test  confirming  Laplace’s  law  for  pressure  difference. 


SOR  method,  the  derivation  method  of  density  with  steep  gradi¬ 
ents,  and  other  refinements,  make  a  stable  and  reliable  simulation 
of  two-phase  flows  with  large  density  differences  possible  [9]. 

To  validate  the  present  LBM  model,  a  liquid  water  droplet  is  ini¬ 
tially  placed  at  the  center  of  a  50  x  50  x  50  (in  lattice  units)  domain 
without  gravity  and  air  current.  According  to  Laplace  law,  when  the 
system  reaches  the  equilibrium  condition  in  the  absence  of  body 
force,  the  pressure  difference  between  interior  and  exterior  of  the 
droplet,  A p,  is  determined  by  the  radius  of  the  droplet  R  and  the 
surface  tension  a  as: 


A  p 


2 

—  Pint  -  Pext  — 


(11) 


To  test  Laplace’s  law,  given  by  Eq.  (11),  the  simulation  was  con¬ 
ducted  for  several  droplet  radii.  The  change  of  pressure  difference 
with  respect  to  2 /R  is  plotted  in  Fig.  2  and  it  exhibits  good  agree¬ 
ment  with  Laplace’s  law.  The  slope  of  the  linear  fit  is  the  interfacial 
tension,  cr,  which  is  about  0.067  N  m-1  for  the  present  droplet  test 
simulation. 

Fig.  3  is  a  schematic  outline  of  the  domain  for  the  cal¬ 
culations  used  in  the  simulations.  The  density  ratio  of  liquid 
to  gas  is  pLl pG  =  847  (pL  =  997  kg  m-3,  pG  =  1.18  kgnrr3),  the  vis¬ 
cosities  of  the  liquid  and  gas  are  fiL  =  8.54  xlO-4 Pas  and 
fiG  =  1 .86  x  1 0-5  Pa  s,  and  the  interfacial  tension  between  water  and 
air  is  cr  =  7.29  x  1 0-2  N  m-1 .  The  time  step  At  is  set  to  2.5  x  1 0-7  s 
and  the  gravitational  acceleration  is  g=0ms-2,  with  other 


0.5  [mm] 


Fig.  3.  Model  for  calculating  the  behavior  of  a  liquid  water  droplet  in  a  gas  channel. 


parameters  zy=  1,  rg  =  1,  Kf=  0.5(Ax)2,  4>L  =  0.092  and  0G  =  0.015. 
The  domain  is  divided  into  a  40  x  20  x  100  cubic  cells  of  0.025  mm 
in  the  x,y,  and  z  directions  with  the  channel  length  2.5  mm,  which  is 
sufficient  for  a  single  droplet  simulation.  The  bottom  of  the  channel 
corresponds  to  the  gas  diffusion  layer  (GDL),  which  is  a  hydropho¬ 
bic  surface  with  an  equilibrium  static  contact  angle  (0S)  of  120°  and 
the  order  parameter  05  is  equal  to  0.045.  The  relationship  between 
the  equilibrium  contact  angle  0S  and  the  order  parameter  05  will  be 
discussed  later.  The  other  three  walls  are  also  hydrophobic  surfaces 
with  static  contact  angles  of  100°  (05  =  0.050). 

The  roughness  of  a  porous  medium  like  GDL  affects  the  droplet 
movement,  but  this  model  assumes  a  smooth  surface  and  so  any 
effect  of  GDL  roughness  is  ignored,  like  in  Ref.  [5].  The  liquid  water 
droplet  is  placed  either  at  a  corner  or  at  the  center  of  the  cell  (Fig.  3 
shows  the  case  with  the  droplet  at  the  corner).  Gravitational  forces 
were  not  considered  in  this  simulation.  A  Poiseuille-like  flow  is 
given  at  the  inlet  of  the  channel,  z  =  0  and  a  free  outflow  condition 
is  used  at  the  outlet  of  the  channel. 

To  evaluate  the  effectiveness  of  the  LBM  model  here  for  wetta¬ 
bility  effects  in  the  simulation,  a  static  droplet  test  was  performed. 
The  static  contact  angle  is  represented  by  the  solid  wall  index  func¬ 
tion  05,  and  in  the  present  contact  angle  simulation,  initially,  a 
liquid  droplet  of  radius  1 0  in  lattice  units,  is  placed  at  the  geometric 
center  of  the  bottom  solid  wall  of  the  40  x  20  x  40  (in  lattice  units) 
computational  domain.  Here,  only  the  bottom  surface  index  func¬ 
tion,  which  is  in  contact  with  the  droplet,  is  changed.  The  value  of  0s 
is  varied  until  the  droplet  reaches  equilibrium  with  an  unchanged 
spherical-cap  shape  achieving  an  unchanged  droplet  shape  with 
different  contact  angles.  Fig.  4  shows  two  contact  angles  obtained 
by  adjusting  the  bottom  wall  index  function  as  well  as  the  density 
contours  of  the  droplet  fluid  at  the  mid-section.  The  obtained  static 
contact  angles  for  different  wall  index  functions  are  plotted  in  Fig.  5, 
showing  the  contact  angle  as  a  function  of  0S.  Fig.  5  shows  that  a 
higher  than  0.055  value  of  0S  gives  rise  to  a  contact  angle  of  less 
than  90°,  indicating  a  hydrophilic  surface.  A  contact  angle  larger 
than  90°  is  formed  when  0S  is  less  than  0.055;  the  0S  =  0.055  is  the 
neutral  wetting  situation. 

Next,  the  dynamic  behavior  of  a  moving  droplet  placed  at  the 
center  of  the  gas  channel  was  also  simulated  in  a  computational 
40  x  40  x  100  cell  domain  with  the  same  lattice  space  interval  as 
in  the  basic  cell  domain  case.  Initially  a  droplet  at  an  equilibrium 
state  is  placed  at  the  geometric  center  of  the  bottom  solid  wall 
and  moves  due  to  the  Poiseuille  gas  flow  with  a  mean  velocity  of 
Uin  =  0.32  ms-1.  Fig.  6(a)  is  an  enlarged  view  of  the  profile  of  the 
moving  droplet  and  the  corresponding  velocity  fields.  The  veloc¬ 
ity  field  shows  the  direction  and  difference  between  the  gas  and 
droplet  velocities  value.  The  deformation  of  the  moving  droplet  and 
its  motion  on  the  hydrophobic  surface  is  clearly  suggested  by  this 
figure.  The  corresponding  mid-section  of  the  absolute  velocity  pro¬ 
file  is  plotted  in  Fig.  6(b).  The  location  ofy/H  =  0.4  corresponds  to 
the  top  of  the  droplet.  The  effect  of  the  flow  velocity  on  the  dynamic 
contact  angle  was  also  studied.  Fig.  7  shows  simulated  droplet 
shapes  at  different  gas  flow  velocities  for  the  basic  case,  where  the 
initial  static  shape  of  the  droplet  is  shown  to  be  a  truncated  sphere 
(Fig.  7  “initial”).  It  clearly  shows  that  droplet  deformation  increases 
as  the  average  gas  flow  velocity  Uin  increases  due  to  the  increases 
in  the  drag  force  acting  on  the  droplet.  For  higher  flow  velocities, 
the  droplet  is  continuously  deformed  and  displays  a  growing  top 
and  spreading  out  as  especially  shown  for  the  case  of  Uin  =  0.8  m  s-1 ; 
this  deformation  will  lead  to  droplet  instability,  and  so  easier  water 
removal  from  the  channel. 

The  droplet  movement  induces  different  advancing  and  reced¬ 
ing  contact  angles,  0A  and  Or.  Fig.  8  presents  the  results  of  the 
LBM  simulations  for  the  contact  angle  hysteresis,  cos0r-cos0a, 
as  a  function  of  the  droplet  velocity  Udr  for  steady  conditions.  The 
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Fig.  4.  Wettability  change  with  different  solid  wall  index  functions:  (a)  0S  =  0.045, 6S  =  120° ;  (b)  0S  =  0.060, 6S  =  70°. 


contact  angle  hysteresis  increases  with  increases  in  Udr.  The  depen¬ 
dence  is  close  to  linear  which  is  in  good  agreement  with  the  work 
of  Hao  et  al.  [16]  based  on  a  macroscopic  force  balance  analyti¬ 
cal  analysis.  The  simulated  linear  functional  dependence  between 
the  contact  angle  hysteresis  and  droplet  terminal  velocity  was  also 
observed  by  the  experiments  by  Kumbur  et  al.  [3].  This  shows  a 
partial  validation  of  the  model  within  the  limited  available  data.  It 
should  be  noted  that  the  LBM  can  estimate  the  relation  between 
the  contact  angle  hysteresis  and  the  droplet  motion  without  any 
experimental  data  about  the  dynamic  contact  angle  of  the  moving 
droplet. 

Fig.  9  shows  the  changes  in  the  movement  of  the  center  of  grav¬ 
ity  of  a  liquid  water  droplet  along  the  channel  and  the  pressure  drop 
in  the  air  flow  in  a  1 .0  mm  wide  channel  for  the  basic  case  (detailed 
in  Section  3.1 ).  The  air  flow  rate  is  equal  to  24  seem  (standard  cubic 
centimeters  per  minute).  This  value  is  very  similar  to  an  actual  fuel 
cell  under  the  following  operation  conditions:  cell  current  density 
7=0.5  A  cm-2,  active  area  of  2  cm2  for  one  100  mm  long  channel, 
and  a  stoichiometric  ratio  of  about  1.4.  The  initial  droplet  position 
is  at  the  center.  The  gas  mean  flow  velocity  Uin  is  0.8  m  s-1 ,  and  the 
time  step  is  2.5  x  10-7  s.  Fig.  9  shows  that  the  velocity  of  the  liquid 
water  droplet,  which  corresponds  to  the  gradient  of  the  moving 
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Fig.  5.  Calculated  static  contact  angle  with  wall  index  function. 


distance,  and  the  pressure  drop  of  air  flow  reach  steady  values.  The 
terminal  droplet  velocity  and  the  pressure  of  the  air  flow  will  be 
used  for  the  evaluation  of  the  draining  performance  in  the  following 
sections. 


Fig.  6.  Velocity  profiles:  (a)  view  of  the  profile  of  the  moving  liquid  droplet  and  the 
corresponding  relative  air  flow  velocities:  (b)  the  corresponding  cross-section  of  the 
absolute  velocity  profile. 
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Fig.  7.  Droplet  behaviors  for  different  mean  air  flow  velocities. 


Fig.  8.  Plot  of  the  contact  angle  hysteresis  vs.  droplet  velocity  relation,  showing  a 
good  agreement  with  Refs.  [3,16]. 


Fig.  9.  Moving  distance  of  a  liquid  water  droplet  and  the  pressure  drop  in  the  gas 
flow. 


3.2.  Evaluation  of  draining  performance 


32  A.  Effect  of  channel  height 

The  simulations  of  droplet  behaviors  for  droplets  of  different 
sizes  placed  at  the  center  of  the  hydrophobic  bottom  wall  were 
performed  for  different  channel  heights  in  channels  with  the  same 
width  and  for  one  gas  flow  rate  corresponding  to  the  basic  condi¬ 
tion  with  a  gas  inlet  velocity  Uin  is  0.8  m  s-1 .  Here,  for  larger  droplet 
volumes,  the  droplet  touches  the  top  wall  for  shallow  channels 
H  =  0.4  mm  and  0.5  mm,  and  for  deeper  channels  H  =  0.7  mm  and 
1.0  mm  the  droplet  touches  the  two  side  walls  of  the  gas  chan¬ 
nel.  In  Fig.  10,  the  cases  where  the  droplets  touch  the  top  wall  are 
marked  with  a  circle  (O)  and  droplets  touching  the  two  side  walls 
are  marked  with  a  dashed  square  (G ).  Fig.  10(a)  shows  the  changes 
in  the  droplet  terminal  velocity,  Udr  vs.  the  droplet  volume  for  chan¬ 
nel  heights  of  0.4-1 .0  mm.  As  the  water  droplet  size  in  the  channel 
increases,  the  droplet  terminal  velocity  also  increases.  When  the 
droplet  touches  the  top  wall  or  the  side  walls,  the  terminal  velocity 
decreases  due  to  the  resistance  arising  from  the  wall.  The  decrease 
in  terminal  velocity  of  the  droplets  will  cause  an  increase  of  air  flow 
pressure.  The  pressure  drop  in  the  gas  channel  is  analyzed  in  terms 
ofpsd  (the  difference  of  actual  pressure  drop  with  a  droplet  present 
and  the  pressure  drop  without  droplets).  The  results  of  the  pres¬ 
sure  drop  in  the  gas  channel  are  plotted  in  Fig.  10(b).  Comparing 
the  cases  where  droplets  either  touch  the  top  wall  or  the  side  walls, 
show  that,  for  the  same  droplet  mass,  deeper  channels  result  in  a 
lower  pressure  drop. 

After  analyzing  the  droplet  velocity  and  the  pressure  drop,  we 
introduced  a  “pumping  efficiency  parameter”,  rjt  defined  as  follows: 


mdr§Udr 

PsdQ 


(12) 


Here,  Q.  is  the  gas  flow  rate  and  mdr  is  the  liquid  droplet’s  mass. 
When  we  consider  the  frictional  work  of  moving  a  droplet  at  a 
velocity  of  Udr,  the  power  is  proportional  to  mdrgUdr.  Thus  the 
pumping  efficiency  has  a  meaning  of  droplet  moving  power  relative 
to  the  pumping  work.  In  other  words,  this  is  related  to  the  parame¬ 
ter  inversely  proportional  to  the  effective  friction  coefficient.  Larger 
pumping  efficiency  indicates  smaller  equivalent  friction  coefficient, 
and  it  corresponds  to  the  better  water  removal  ability  for  the  same 
compressor  work.  In  general,  a  low  pressure  difference  across  the 
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Fig.  11.  Changes  in  (a)  droplet  terminal  velocities  and  (b)  pumping  efficiency  at 
corners  for  different  droplet  masses  with  different  channel  heights. 

to  the  case  where  the  droplet  touches  the  corner  or  the  top  wall. 
Considering  these  results,  the  channel  height  of  0.5  mm  may  be  con¬ 
cluded  to  be  superior  for  draining  the  droplets  since  the  drainage 
speed  is  also  high.  This  optimum  channel  shape  was  also  reported 
in  the  experimental  work  of  Akhtar  et  al.  from  a  different  viewpoint 
[17].  They  concluded  from  the  experiments  of  a  droplet  detachment 
in  a  channel  that  the  rectangular  shaped  channel  with  a  width  of 
1  mm  and  a  depth  of  0.5  mm  is  found  to  exhibit  best  water  removal 
properties  at  a  reasonable  pressure  drop. 


Fig.  10.  Changes  in  (a)  droplet  terminal  velocities,  (b)  pressure  drop  and  (c)  pumping 
efficiency  for  different  droplet  masses  with  different  channel  heights. 

flow  field  is  desired  because  of  lower  auxiliary  energy  demand,  e.g. 
for  air  compression.  The  pumping  efficiency  is  plotted  vs.  droplet 
size  in  Fig.  10(c).  The  plots  show  two  regions:  the  first  corresponds 
to  the  case  when  the  droplet  does  not  touch  side  walls  or  top  wall 
and  the  second  region  is  for  the  case  where  the  droplet  touches 
the  top  wall  or  the  two  side  walls  of  the  gas  channel.  When  the 
droplet  does  not  touch  the  walls,  the  pumping  efficiency  increases 
with  increases  in  channel  height  leading  to  a  slowing  of  the  speed 
of  the  droplet  since  the  flow  rate  of  the  gas  is  constant.  For  the  sec¬ 
ond  region,  where  the  droplet  touches  the  top  wall,  for  H  =  0.4  mm 
and  H  =  0.5mm,  the  pumping  efficiency  is  only  little  affected  by 
the  touching  of  the  top  wall  effect.  However,  when  the  droplet 
touches  the  side  walls  in  the  deeper  gas  channels,  H  =  0.7mm 
and  H=1.0mm,  the  pumping  efficiency  is  dramatically  decreased 
because  of  the  higher  resistance  exerted  by  the  two  side  walls  on 
the  droplet  motion  and  the  pumping  efficiency  becomes  similar 


3.2.2.  Effect  of  droplet  position 

The  initial  droplet  position  is  also  an  important  factor  in  the 
pump  work,  and  the  simulation  of  the  water  behavior,  where  the 
droplet  is  placed  at  the  corner  of  the  channel  is  shown  in  Fig.  11(a), 
and  shows  remarkable  differences  from  the  droplet  velocity 
values  of  the  case  when  droplet  is  placed  at  the  center  (Fig.  10(a)). 
When  the  droplet  is  placed  at  the  corner,  the  water  droplet  is  sub¬ 
ject  to  a  high  resistance  from  the  contact  walls  and  it  is  far  from 
the  Poiseuille  flow  mainstream  (center  of  the  channel  where  the 
velocity  is  the  highest).  In  this  situation,  the  droplet  velocity  is  much 
slower,  which  explains  the  large  differences  between  Fig.  10(a)  and 
Fig.  11(a).  For  larger  droplet  sizes  and  when  the  droplet  is  placed 
at  corner  and  touching  the  top  wall,  the  terminal  droplet  velocity 
is  slightly  increased  due  to  the  larger  volume  and  as  more  of  the 
droplet  volume  is  nearer  the  gas  channel  center  and  mainstream 
flow,  which  results  in  a  larger  pressure  drop  and  the  pumping 
efficiency  increases  slightly  as  shown  in  Fig.  11(b).  Increasing  the 
droplet  mass  further  changes  the  droplet  shape  to  a  liquid  water 
film. 
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Fig.  12.  Changes  in  (a)  droplet  terminal  velocities  and  (b)  pumping  efficiency  for 
hydrophobic  and  hydrophilic  channel  walls  for  a  droplet  at  a  corner  H  =  0.4  mm  and 
1.0  mm  for  different  droplet  masses. 


3.2.3.  Effect  of  wettability 

In  this  study,  the  wettability  effect  is  considered  in  the  index 
function  of  the  solid  wall.  The  dynamic  simulation  for  both 
hydrophobic  and  hydrophilic  separators  (static  contact  angle 
60°,  0s  =  0.060)  and  the  bottom  surface,  which  corresponds 
to  the  GDL  hydrophobic  as  in  previous  simulations,  was  con¬ 
ducted  for  the  case  when  the  droplet  is  placed  at  the  corner 
of  the  gas  channel.  Fig.  12(a)  displays  the  droplet  terminal 
velocity  change  with  droplet  mass  for  deeper  and  shallower  chan¬ 
nels,  H=1.0mm  and  H  =  0.4mm,  for  both  the  hydrophobic  and 
hydrophilic  cases  when  the  droplet  is  placed  at  a  corner.  For  both 
the  deeper  and  shallower  channels,  the  droplet  velocity  is  rela¬ 
tively  high  for  a  hydrophobic  channel  walls  than  for  a  hydrophilic 
walls.  The  corresponding  pumping  efficiencies  are  plotted  in 
Fig.  12(b).  It  shows  that  in  the  shallow  gas  channel,  the  pumping 
efficiency  is  less  affected  by  the  wettability  and  the  droplet  is  able 
to  maintain  a  relatively  high  velocity  that  will  result  in  a  high 
“drainage  speed”.  In  both  cases  a  hydrophobic  separator  gives  the 
larger  pumping  efficiency. 

Fig.  13  shows  a  three-dimensional  views  of  droplet  behav¬ 
ior  placed  at  corner  of  an  initial  radius  of  16  lattice  units  and  a 
channel  height  of  0.4  mm  for  both  hydrophobic  and  hydrophilic 
channel  walls.  In  this  case  the  droplet  touches  the  top  wall. 
Fig.  13  shows  that  for  the  hydrophilic  separator  case  the  droplet 
attaches  to  the  top  wall  and  the  contact  area  with  GDL  is 
smaller  than  with  the  hydrophobic  separator.  This  is  very  advan¬ 
tageous  for  the  fuel  cell  performance.  A  minimum  bottom  liquid 


Fig.  13.  A  three-dimensional  view  of  a  droplet  (radius  14  lattice  units)  H  =  0.4mm 
(a)  with  hydrophobic  channel  walls  and  (b)  hydrophilic  channel  walls. 


contact  area  makes  more  space  available  for  the  gas  diffu¬ 
sion. 

3.3.  Design  concept 

Understanding  the  basic  phenomena  and  dynamic  behavior  of  a 
liquid  water  droplet  in  a  single  gas  channel  here  should  be  extended 
to  more  complex  and  larger  scale  fuel  cell  simulation  like  Ref.  [18]. 
This  section  presents  some  discussion  about  the  application  of  the 
simulation  results  for  design  concept  of  a  fuel  cell  channel.  In  the 
large  scale  channel,  a  drained  moving  droplet  unites  with  other 
droplets,  grows  larger  and  fills  with  the  channel.  The  plugging 
induces  a  drastic  increase  in  pressure  drop  of  the  air  flow,  but  the 
droplet  is  drained  immediately  with  similar  velocity  to  the  air  flow. 
Therefore,  it  is  important  that  the  droplet  is  moved  before  the  plug¬ 
ging  with  high  velocity,  high  pumping  efficiency  and  small  contact 
area  with  the  GDL.  The  simulation  results  showed  that  the  droplet 
velocity  with  shallower  channels  is  kept  higher  and  the  pumping 
efficiency  becomes  less  dependent  on  the  droplet  locations  with 
shallower  channels  as  the  droplet  volume  increases.  This  may  lead 
to  the  same  conclusion  in  large  scale  channels  that  shallower  chan¬ 
nels  with  about  0.5  mm  height  are  superior  to  deeper  channels.  It 
was  also  shown  that  a  hydrophilic  channel  wall,  relative  to  the  GDL 
wall,  is  better  for  minimum  liquid  contact  area  with  the  GDL  for 
larger  droplet  volume. 

4.  Conclusions 

The  paper  investigates  water  droplet  behavior  in  fuel  cell  chan¬ 
nels  with  the  same  width  at  the  same  gas  flow  rate  condition,  but 
for  different  gas  channel  heights,  droplet  positions  and  gas  channel 
wall  wettability,  using  LBM  simulation.  The  water  drainage  per¬ 
formance  can  be  characterized  by  the  two  parameters  water  flow 
rate  and  pump  work,  which  are  expressed  by  the  droplet  velocity 
and  pumping  efficiency  in  the  paper.  The  results  obtained  in  the 
analysis  may  be  summarized  as  follows: 

1 .  Droplet  velocity  significantly  decreases  when  the  droplets  touch 
the  corners  or  the  top  wall  compared  to  the  case  where  the 
droplet  locates  on  the  center  of  the  GDL  surface  without  touching 
the  side  or  top  walls. 

2.  Deeper  channels  give  better  drain  efficiency  than  shallower 
channels,  but  the  efficiency  differences  become  small  when  the 
droplet  touches  to  the  corner  or  the  top  wall.  As  the  droplet 
velocity,  i.e.  the  draining  flow  rate,  becomes  higher  and  the 
pumping  efficiency  becomes  less  dependent  on  the  droplet 
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locations  with  shallower  channels,  shallower  channels  are  supe¬ 
rior  to  deeper  channels  when  the  pump  work  involved  are 
similar. 

3.  Compared  to  hydrophobic  channel  walls,  hydrophilic  walls  may 
result  in  better  gas  transport  characteristics,  as  the  liquid  water 
is  drawn  up  on  the  channel  wall  from  the  GDL  surface  to 
leave  more  open  area  available  for  gas  transport  to  the  GDL 
Hydrophilic  walls  result  in  a  larger  pressure  drop  and  lower 
draining  flow  rates  than  hydrophobic  walls,  however,  the  dif¬ 
ferences  are  smaller  with  shallower  channels. 
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